// Copyright (C) 2011  Andrew H. Chan
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
// 
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.

#ifndef _NTSEQ_H__
#define _NTSEQ_H__

#include <iostream>
#include <vector>
#include <assert.h>

typedef unsigned char uchar;
typedef unsigned long long uint64;

// no ability to remove bases
// no ability to insert bases, except at end of sequence
class ntSeq {
    public:
        static char invBM[4];
        static uchar BM[52];

        ntSeq() : numBases(0), bases(0) {}
        
        ntSeq(std::string seq) : numBases(0), bases(0) {
            for(uint64 ii=0; ii<seq.size(); ++ii) {
                append(bm(seq[ii]));
            }
        }

        inline uchar bm(char base) { return BM[base-65]; }

        inline void append(uchar base) {
            if(!(numBases%32)) bases.push_back(0);
            set(numBases, base);
            ++numBases;
        }

        inline void set(size_t index, uchar base) {
            assert(index<size()+1);
            assert(base>=0&&base<4);
            size_t majorIndex = index/32;
            uint64 minorShift = 2*(31 - index%32); // bases are stored starting from MSB
            uint64 mask = (uint64)0x3 << minorShift;
            uint64 baseMask = (uint64)base << minorShift;
            bases[majorIndex] &= ~mask; // zero out position of interest
            bases[majorIndex] |= baseMask; // set position of interest
        }

        inline uchar get(size_t index) const {
            assert(index<size());
            size_t majorIndex = index/32;
            uint64 minorShift = 2*(31 - index%32); // bases are stored starting from MSB
            uint64 mask = (uint64)0x3 << minorShift;
            return (uchar)((bases[majorIndex] & mask) >> minorShift);
        }

        void resize(size_t newSize) {
            if(newSize==0) {
                bases.resize(0);
            } else {
                bases.resize((newSize-1)/32+1);
            }
            numBases = newSize;
        }

        inline size_t size() const {
            return numBases;
        }

        std::string getSuffix(uint64 pos) const {
            // Should be able to handle pos == size() (empty suffix)
            assert(pos<=size());
            std::string tmp;
            for(uint64 ii=pos; ii<size(); ++ii) tmp += invBM[get(ii)];
            return tmp;
        }

        size_t numBases;
        std::vector<uint64> bases;
};

#endif
